#pragma once
#include "gtest/gtest.h"
#include "DataStructure/Common/Real.h"
#include "DBManager/DBServiceFactory.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Coordinate/ICORD2RService.h"
#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "../../Factory/ElementCalculatorFactory.h"
#include "../BarCalculator.h"

using SG::Algebra::Matrixd;
using SG::Algebra::Vector3D;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::Common::Real;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::CoordinateProcessor;
using SG::FEMSolver::Element::ElementCalculatorFactory;
using SG::FEMSolver::Element::BarCalculator;

// ElementCalculatorFactory 是否正确创建 Bar 单元计算类
TEST (ElementCalculatorFactoryTest, CreateBarCalculator)
{
    ElementCalculatorFactory::Register (ElementTypes::ELETYPE_Bar, BarCalculator::Create);

    // Test data
    const auto calculatorPtr = ElementCalculatorFactory::Create (ElementTypes::ELETYPE_Bar);
    const auto barCalculatorPtr = dynamic_cast<BarCalculator*> (calculatorPtr.get ());

    // Check
    ASSERT_TRUE (nullptr != barCalculatorPtr);
}

extern std::shared_ptr<DBServiceFactory> DB_SERVICEFACTORY;
class BarCalculatorTest : public ::testing::Test
{
protected:
    static void SetUpTestSuite ()
    {
        DB_SERVICEFACTORY->ClearDatabase ();

        // 直角坐标系1
        CORD2RData rectangularCoordSys{1, 0, {1.0, 2.0, 3.0}, {2.0, 2.0, 3.0}, {1.0, 3.0, 3.0}};
        DB_SERVICEFACTORY->get<ICORD2RService> ()->Add (rectangularCoordSys);

        // 直角坐标系2
        rectangularCoordSys = {2, 0, {4.0, 2.0, 1.0}, {4.0, 3.0, 1.0}, {4.0, 2.0, 2.0}};
        DB_SERVICEFACTORY->get<ICORD2RService> ()->Add (rectangularCoordSys);

        // 生成节点排序
        auto& assembly = AssemblyEigen::GetInstance ();
        assembly.Clear ();
        Subcase subcase;
        assembly.Initialize (DB_SERVICEFACTORY, subcase, {1});

        // 计算坐标系原点和基矢量
        auto& coordProcessor = CoordinateProcessor::GetInstance ();
        coordProcessor.Clear ();
        coordProcessor.Initialize (DB_SERVICEFACTORY);
    }

    void SetUp () override
    {
        // 小量放大系数
        factor = 1.0E12;
    }

    Real factor;
};

TEST_F (BarCalculatorTest, ComputeOffset)
{
    // Common Input data
    const Matrixd TO (3, 3, {
        0.23076923076923078,  0.30769230769230771,  0.92307692307692313,
        0.96071401390504119, -0.22238750321875955, -0.16604933573667377,
        0.15418866889833996,  0.92513201339003959, -0.34692450502126487
    });

    // Case 1: OFFT 为 BOG，节点位移坐标系 CD = 0
    {
        // Input data
        const std::vector<BeamCoordOptions> offsetFlagSet{BeamCoordOptions::OPTION_B, BeamCoordOptions::OPTION_O, BeamCoordOptions::OPTION_G};
        const IdArray dispCoordIdSet{0, 0};
        const std::vector<Matrixd> TGSet (2);
        std::vector<Vector3D> offsetSetTest{{2.0, 8.0, 4.0}, {6.0, 9.0, 11.0}};

        // Test data
        std::vector<Vector3D> offsetSetGTest (2);
        BarCalculator::ComputeOffset (offsetFlagSet, TO, dispCoordIdSet, TGSet, _OUT offsetSetTest, _OUT offsetSetGTest);

        // Target data
        std::vector<Vector3D> offsetSetTarget{{8.7640052483721504, 2.5368126431946973, -0.86993885982460339}, {6.0, 9.0, 11.0}};
        std::vector<Vector3D> offsetSetGTarget = offsetSetTarget;
        
        // Check
        ASSERT_TRUE (offsetSetTarget == offsetSetTest);
        ASSERT_TRUE (offsetSetGTarget == offsetSetGTest);
    }

    // Case 2：OFFT 为 BGO，节点位移坐标系 CD > 0
    {
        // Input data
        const std::vector<BeamCoordOptions> offsetFlagSet{BeamCoordOptions::OPTION_B, BeamCoordOptions::OPTION_G, BeamCoordOptions::OPTION_O};
        const IdArray dispCoordIdSet{2, 1};
        const std::vector<Matrixd> TGSet{
            Matrixd (3, 3, {
                0.0, 0.0, 1.0,
                1.0, 0.0, 0.0,
                0.0, 1.0, 0.0
            }),
            Matrixd (3, 3, {
                0.0, 1.0, 0.0,
                0.0, 0.0, 1.0,
                1.0, 0.0, 0.0
            })
        };
        std::vector<Vector3D> offsetSetTest{{2.0, 8.0, 4.0}, {6.0, 9.0, 11.0}};

        // Test data
        std::vector<Vector3D> offsetSetGTest (2);
        BarCalculator::ComputeOffset (offsetFlagSet, TO, dispCoordIdSet, TGSet, _OUT offsetSetTest, _OUT offsetSetGTest);

        // Target data
        std::vector<Vector3D> offsetSetTarget{{8.0, 4.0, 2.0}, {11.727116867642495, 10.021118464475446, 0.22784796159756127}};
        std::vector<Vector3D> offsetSetGTarget{{2.0, 8.0, 4.0}, {10.021118464475446, 0.22784796159756127, 11.727116867642495}};

        // Check
        ASSERT_TRUE (offsetSetTarget == offsetSetTest);
        ASSERT_TRUE (offsetSetGTarget == offsetSetGTest);
    }
}

TEST_F (BarCalculatorTest, ComputeElementTransWithoutOffset)
{
    // Input data
    std::vector<NodeData> nodeSet (2);
    NodeData orienNode;

    CBARData element;
    element.m_id   = 1;
    element.m_g[0] = 1;
    element.m_g[1] = 2;
    element.m_go   = 3;

    // Target data
    const Real LTarget{sqrt (2.0)};
    const Matrixd TEleTarget (12, 12, {
    0.0, 0.70710678,  0.70710678, 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    1.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.70710678, -0.70710678, 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.70710678,  0.70710678,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 1.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.70710678, -0.70710678,  0.0       , 0.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.70710678, 0.0, 0.70710678,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.0       , 1.0, 0.0       ,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , -0.70710678, 0.0, 0.70710678,  0.0       , 0.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.70710678, 0.0, 0.70710678,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       ,  0.0       , 1.0, 0.0       ,
    0.0, 0.0       ,  0.0       , 0.0, 0.0       ,  0.0       ,  0.0       , 0.0, 0.0       , -0.70710678, 0.0, 0.70710678
    });

    Real LTest;
    Matrixd T (3, 3);
    Matrixd TEleTest (12, 12);

    // G1坐标系为basic（id = 0）；G0坐标系为basic
    {
        // Input data
        nodeSet = {{1, 0, {0.0, 0.0, 0.0}, 2}, {2, 2, {-1.0, -3.0, -1.0}, 1}};
        orienNode = {3, 0, {0.0, 0.0, 2.0}, 1};

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }

    // G1坐标系为local（id > 0）；G0坐标系为local
    {
        // Input data
        nodeSet = {{1, 1, {-2.0, -3.0, -1.0}, 2}, {2, 2, {-1.0, -3.0, -1.0}, 1}};
        orienNode = {3, 2, {1.0, -4.0, -2.0}, 1};

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }
}

TEST_F (BarCalculatorTest, ComputeElementTransWithOffset)
{
    // Input data
    const std::vector<NodeData> nodeSet{{1, 1, {-2.0, -3.0, -1.0}, 2}, {2, 2, {-1.0, -3.0, -1.0}, 1}};
    NodeData orienNode;

    CBARData element;
    element.m_id   = 1;
    element.m_g[0] = 1;
    element.m_g[1] = 2;
    element.m_w1b  = 0.0;
    element.m_w2b  = 1.0;
    element.m_w3b  = 0.0;

    // Target data
    const Real LTarget{sqrt (2.0)};
    const Matrixd TEleTarget (12, 12, {
    0.70710678, 0.0,  0.70710678,  0.70710678, 0.0, -0.70710678,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.70710678, 0.0, -0.70710678, -0.70710678, 0.0, -0.70710678,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.0       , 1.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.70710678, 0.0,  0.70710678,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.70710678, 0.0, -0.70710678,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.0       , 1.0,  0.0       ,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.70710678, 0.70710678, 0.0,  0.0       , 0.0       , -0.70710678,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       , -0.70710678, 0.70710678, 0.0,  0.0       , 0.0       ,  0.70710678,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.0       , 0.0       , 1.0,  1.0       , 0.0       ,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.0       , 0.0       , 0.0,  0.70710678, 0.70710678,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.0       , 0.0       , 0.0, -0.70710678, 0.70710678,  0.0       ,
    0.0       , 0.0,  0.0       ,  0.0       , 0.0,  0.0       ,  0.0       , 0.0       , 0.0,  0.0       , 0.0       ,  1.0
    });

    Real LTest;
    Matrixd T (3, 3);
    Matrixd TEleTest (12, 12);

    // 朝向矢量, GOO
    {
        // Input data
        element.m_flag = 21;  // GOO
        element.m_x1   = 2.0;
        element.m_x2   = 0.0;
        element.m_x3   = 0.0;
        element.m_go   = 0;
        element.m_w1a  = 0.5 * sqrt (2.0);
        element.m_w2a  = 0.0;
        element.m_w3a  = 0.5 * sqrt (2.0);

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }

    // 朝向矢量, BGO
    {
        // Input data
        element.m_flag = 81;  // BGO
        element.m_x1   = 0.0;
        element.m_x2   = 0.0;
        element.m_x3   = 2.0;
        element.m_w1a  = 0.0;
        element.m_w2a  = 1.0;
        element.m_w3a  = 0.0;

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }

    // 朝向节点G0，坐标系为local（id > 0）；偏置矢量GO
    {
        // Input data
        orienNode = {3, 2, {1.0, -4.0, -2.0}, 1};
        element.m_go = 3;

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }

    // 朝向节点G0，坐标系为basic（id = 0）；偏置矢量OO
    {
        // Input data
        orienNode = {3, 0, {0.0, 0.0, 2.0}, 1};
        element.m_flag = 21;  // GOO
        element.m_w1a  = 0.5 * sqrt (2.0);
        element.m_w2a  = 0.0;
        element.m_w3a  = 0.5 * sqrt (2.0);

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT LTest, _OUT T, _OUT TEleTest);

        // Check
        ASSERT_TRUE (LTarget == LTest);
        ASSERT_TRUE (TEleTarget == TEleTest);
    }
}

TEST_F (BarCalculatorTest, ComputeLocalStiffnessWithoutShear)
{
    // Input data
    const Real L{10.0};

    IsotropicMaterial material;
    material.m_E  = 1.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2 * (1 + material.m_nu));

    // Section constant
    PBARData property;
    property.m_A   = 6.0;
    property.m_I1  = 2.0;
    property.m_I2  = 4.5;
    property.m_I12 = 0.0;
    property.m_J   = 4.5 + 2.0;
    property.m_K1  = 0.0;
    property.m_K2  = 0.0;

    // Test data
    Matrixd kgTest (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, kgTest);

    // Target data
    const Matrixd kgTarget (12, 12, {
     0.60000000,  0.00000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000, -0.60000000,  0.00000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.24000000E-01,  0.00000000    ,  0.00000000,  0.00000000,  0.12000000,  0.00000000, -0.24000000E-01,  0.00000000    ,  0.00000000,  0.00000000,  0.12000000,
     0.00000000,  0.00000000    ,  0.54000000E-01,  0.00000000, -0.27000000,  0.00000000,  0.00000000,  0.00000000    , -0.54000000E-01,  0.00000000, -0.27000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000    ,  0.25000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000    , -0.25000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    , -0.27000000    ,  0.00000000,  1.80000000,  0.00000000,  0.00000000,  0.00000000    ,  0.27000000    ,  0.00000000,  0.90000000,  0.00000000,
     0.00000000,  0.12000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.80000000,  0.00000000, -0.12000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.40000000,
    -0.60000000,  0.00000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.60000000,  0.00000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,
     0.00000000, -0.24000000E-01,  0.00000000    ,  0.00000000,  0.00000000, -0.12000000,  0.00000000,  0.24000000E-01,  0.00000000    ,  0.00000000,  0.00000000, -0.12000000,
     0.00000000,  0.00000000    , -0.54000000E-01,  0.00000000,  0.27000000,  0.00000000,  0.00000000,  0.00000000    ,  0.54000000E-01,  0.00000000,  0.27000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000    , -0.25000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000    ,  0.25000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    , -0.27000000    ,  0.00000000,  0.90000000,  0.00000000,  0.00000000,  0.00000000    ,  0.27000000    ,  0.00000000,  1.80000000,  0.00000000,
     0.00000000,  0.12000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.40000000,  0.00000000, -0.12000000    ,  0.00000000    ,  0.00000000,  0.00000000,  0.80000000
    });

    // Check
    ASSERT_TRUE (kgTarget == kgTest);
}

TEST_F (BarCalculatorTest, ComputeLocalStiffnessWithI2eq0I12eq0)
{
    // Input data
    const Real L{10.0};

    IsotropicMaterial material;
    material.m_E  = 1.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2 * (1 + material.m_nu));

    // Section constant
    PBARData property;
    property.m_A   = 6.0;
    property.m_I1  = 2.0;
    property.m_I2  = 0.0;
    property.m_I12 = 0.0;
    property.m_J   = 4.5 + 2.0;
    property.m_K1  = 0.0;
    property.m_K2  = 0.0;

    // Test data
    Matrixd kgTest (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, kgTest);

    // Target data
    const Matrixd kgTarget (12, 12, {
     0.60000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000, -0.60000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.24000000E-01,  0.00000000,  0.00000000,  0.00000000,  0.12000000,  0.00000000, -0.24000000E-01,  0.00000000,  0.00000000,  0.00000000,  0.12000000,
     0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000,  0.25000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000, -0.25000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.12000000    ,  0.00000000,  0.00000000,  0.00000000,  0.80000000,  0.00000000, -0.12000000    ,  0.00000000,  0.00000000,  0.00000000,  0.40000000,
    -0.60000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.60000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000, -0.24000000E-01,  0.00000000,  0.00000000,  0.00000000, -0.12000000,  0.00000000,  0.24000000E-01,  0.00000000,  0.00000000,  0.00000000, -0.12000000,
     0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000, -0.25000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000,  0.25000000,  0.00000000,  0.00000000,
     0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.00000000    ,  0.00000000,  0.00000000,  0.00000000,  0.00000000,
     0.00000000,  0.12000000    ,  0.00000000,  0.00000000,  0.00000000,  0.40000000,  0.00000000, -0.12000000    ,  0.00000000,  0.00000000,  0.00000000,  0.80000000
    });

    // Check
    ASSERT_TRUE (kgTarget == kgTest);
}

TEST_F (BarCalculatorTest, ComputeGlobalStiffnessWithIyz)
{
    // Input data
    const std::vector<NodeData> nodeSet{{7, 0, {0.0, 0.0, 0.0}, 0}, {8, 0, {100.0, 0.0, 0.0}, 0}};
    const NodeData orienNode;

    CBARData element;
    element.m_id   = 1;
    element.m_pId  = 1;
    element.m_g[0] = 7;
    element.m_g[1] = 8;
    element.m_flag = 1;  // GGG
    element.m_x1   = 0.176992;
    element.m_x2   = 3.239967;
    element.m_x3   = 1.86314;

    IsotropicMaterial material;
    material.m_id = 1;
    material.m_E  = 210000.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2.0 * (1.0 + material.m_nu));

    PBARData property;
    property.m_id  = 1;
    property.m_mId = 1;

    Real L;
    Matrixd T (3, 3);
    Matrixd TEle (12, 12);
    Matrixd kgTest (12, 12);
    Matrixd kgTarget (12, 12);

    // 对称截面Iyz = 0.0
    {
        // Input data
        property.m_A   = 7.0;
        property.m_I1  = 7.145833;
        property.m_I2  = 2.333333;
        property.m_I12 = 0.0;
        property.m_J   = 6.003188;
        property.m_K1  = 0.8333333;
        property.m_K2  = 0.8333333;

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
        BarCalculator::ComputeStiffness (L, material, property, _OUT kgTest);
        kgTest = TEle.Transpose () * kgTest * TEle;

        // Target data
        kgTarget = Matrixd (12, 12, {
         1.47000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.47000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00,  1.49404050E+01,  5.21439250E+00,  0.00000000E+00, -2.60719625E+02,  7.47020249E+02,  0.00000000E+00, -1.49404050E+01, -5.21439250E+00,  0.00000000E+00, -2.60719625E+02,  7.47020249E+02,
         0.00000000E+00,  5.21439250E+00,  8.87120161E+00,  0.00000000E+00, -4.43560080E+02,  2.60719625E+02,  0.00000000E+00, -5.21439250E+00, -8.87120161E+00,  0.00000000E+00, -4.43560080E+02,  2.60719625E+02,
         0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.84872877E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.84872877E+03,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -2.60719625E+02, -4.43560080E+02,  0.00000000E+00,  2.95894603E+04, -1.74033599E+04,  0.00000000E+00,  2.60719625E+02,  4.43560080E+02,  0.00000000E+00,  1.47665477E+04, -8.66860262E+03,
         0.00000000E+00,  7.47020249E+02,  2.60719625E+02,  0.00000000E+00, -1.74033599E+04,  4.98458047E+04,  0.00000000E+00, -7.47020249E+02, -2.60719625E+02,  0.00000000E+00, -8.66860262E+03,  2.48562202E+04,
        -1.47000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.47000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -1.49404050E+01, -5.21439250E+00,  0.00000000E+00,  2.60719625E+02, -7.47020249E+02,  0.00000000E+00,  1.49404050E+01,  5.21439250E+00,  0.00000000E+00,  2.60719625E+02, -7.47020249E+02,
         0.00000000E+00, -5.21439250E+00, -8.87120161E+00,  0.00000000E+00,  4.43560080E+02, -2.60719625E+02,  0.00000000E+00,  5.21439250E+00,  8.87120161E+00,  0.00000000E+00,  4.43560080E+02, -2.60719625E+02,
         0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.84872877E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.84872877E+03,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -2.60719625E+02, -4.43560080E+02,  0.00000000E+00,  1.47665477E+04, -8.66860262E+03,  0.00000000E+00,  2.60719625E+02,  4.43560080E+02,  0.00000000E+00,  2.95894603E+04, -1.74033599E+04,
         0.00000000E+00,  7.47020249E+02,  2.60719625E+02,  0.00000000E+00, -8.66860262E+03,  2.48562202E+04,  0.00000000E+00, -7.47020249E+02, -2.60719625E+02,  0.00000000E+00, -1.74033599E+04,  4.98458047E+04
         });

        // Check
        ASSERT_TRUE (kgTarget == kgTest);
    }

    // 非对称截面Iyz != 0.0
    {
        // Input data
        property.m_A   = 11.0;
        property.m_I1  = 28.95644;
        property.m_I2  = 8.303030;
        property.m_I12 = -7.159091;
        property.m_J   = 14.60417;
        property.m_K1  = 0.3181818;
        property.m_K2  = 0.4545455;

        // Test data
        BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
        BarCalculator::ComputeStiffness (L, material, property, _OUT kgTest);
        kgTest = TEle.Transpose () * kgTest * TEle;

        // Target data
        kgTarget = Matrixd (12, 12, {
         2.31000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.31000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00,  7.56289955E+01,  1.34173593E+01,  0.00000000E+00, -6.70867967E+02,  3.78144977E+03,  0.00000000E+00, -7.56289955E+01, -1.34173593E+01,  0.00000000E+00, -6.70867967E+02,  3.78144977E+03,
         0.00000000E+00,  1.34173593E+01,  1.82648689E+01,  0.00000000E+00, -9.13243447E+02,  6.70867967E+02,  0.00000000E+00, -1.34173593E+01, -1.82648689E+01,  0.00000000E+00, -9.13243447E+02,  6.70867967E+02,
         0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.17956758E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.17956758E+04,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -6.70867967E+02, -9.13243447E+02,  0.00000000E+00,  6.08828965E+04, -4.47245311E+04,  0.00000000E+00,  6.70867967E+02,  9.13243447E+02,  0.00000000E+00,  3.04414482E+04, -2.23622656E+04,
         0.00000000E+00,  3.78144977E+03,  6.70867967E+02,  0.00000000E+00, -4.47245311E+04,  2.52096652E+05,  0.00000000E+00, -3.78144977E+03, -6.70867967E+02,  0.00000000E+00, -2.23622656E+04,  1.26048326E+05,
        -2.31000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.31000000E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -7.56289955E+01, -1.34173593E+01,  0.00000000E+00,  6.70867967E+02, -3.78144977E+03,  0.00000000E+00,  7.56289955E+01,  1.34173593E+01,  0.00000000E+00,  6.70867967E+02, -3.78144977E+03,
         0.00000000E+00, -1.34173593E+01, -1.82648689E+01,  0.00000000E+00,  9.13243447E+02, -6.70867967E+02,  0.00000000E+00,  1.34173593E+01,  1.82648689E+01,  0.00000000E+00,  9.13243447E+02, -6.70867967E+02,
         0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.17956758E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.17956758E+04,  0.00000000E+00,  0.00000000E+00,
         0.00000000E+00, -6.70867967E+02, -9.13243447E+02,  0.00000000E+00,  3.04414482E+04, -2.23622656E+04,  0.00000000E+00,  6.70867967E+02,  9.13243447E+02,  0.00000000E+00,  6.08828965E+04, -4.47245311E+04,
         0.00000000E+00,  3.78144977E+03,  6.70867967E+02,  0.00000000E+00, -2.23622656E+04,  1.26048326E+05,  0.00000000E+00, -3.78144977E+03, -6.70867967E+02,  0.00000000E+00, -4.47245311E+04,  2.52096652E+05
         });

        // Check
        ASSERT_TRUE (kgTarget == kgTest);
    }
}

// 刚度矩阵的自由度释放
TEST_F (BarCalculatorTest, ComputeLocalStiffnessWithPin)
{
    // Input data
    CBARData element;
    element.m_pb = 1;

    const Real L{100.0};

    IsotropicMaterial material;
    material.m_E  = 210000.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2.0 * (1.0 + material.m_nu));

    PBARData property;
    property.m_A   = 7.0;
    property.m_I1  = 7.145833;
    property.m_I2  = 2.333333;
    property.m_I12 = 0.0;
    property.m_J   = 6.003188;
    property.m_K1  = 0.833333;
    property.m_K2  = 0.833333;

    // Test data
    Matrixd kgTest (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kgTest);
    std::vector<std::size_t> freeDofSet;
    for (const auto& dof : element.GetDofSetOfPaRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof));
    }
    for (const auto& dof: element.GetDofSetOfPbRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof) + 6);
    }
    BarCalculator::ReleaseStiffnessDof (freeDofSet, _OUT kgTest);

    // Target data
    const Matrixd kgTarget (12, 12, {
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  1.79389365E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.96946826E+02,  0.00000000E+00, -1.79389365E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.96946826E+02,
    0.00000000E+00,  0.00000000E+00,  5.87267007E+00,  0.00000000E+00, -2.93633503E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -5.87267007E+00,  0.00000000E+00, -2.93633503E+02,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.84872877E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.84872877E+03,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.93633503E+02,  0.00000000E+00,  1.95816745E+04,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.93633503E+02,  0.00000000E+00,  9.78167586E+03,  0.00000000E+00,
    0.00000000E+00,  8.96946826E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.98535906E+04,  0.00000000E+00, -8.96946826E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.98410920E+04,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00, -1.79389365E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -8.96946826E+02,  0.00000000E+00,  1.79389365E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -8.96946826E+02,
    0.00000000E+00,  0.00000000E+00, -5.87267007E+00,  0.00000000E+00,  2.93633503E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.87267007E+00,  0.00000000E+00,  2.93633503E+02,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.84872877E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.84872877E+03,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.93633503E+02,  0.00000000E+00,  9.78167586E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.93633503E+02,  0.00000000E+00,  1.95816745E+04,  0.00000000E+00,
    0.00000000E+00,  8.96946826E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.98410920E+04,  0.00000000E+00, -8.96946826E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.98535906E+04
    });

    // Check
    ASSERT_TRUE (kgTarget == kgTest);
}

// 单元坐标系下的单元结构质量矩阵
TEST_F (BarCalculatorTest, ComputeLocalMass)
{
    // Input data
    const Real L{10.0};

    IsotropicMaterial material;
    material.m_rho = 7.85E-9;

    PBARData property;
    property.m_A   = 7.0;
    property.m_I1  = 7.145833;
    property.m_I2  = 2.333333;
    bool hasTorsion{false};
    bool isBARMASS{false};

    // Test data
    auto lumpedMTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    hasTorsion = true;
    auto lumpedMTorsionTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    auto coupledMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    isBARMASS = true;
    auto consistentMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);

    // Target data
    Matrixd lumpedMTarget (12, 12);
    lumpedMTarget (0, 0) = lumpedMTarget (1, 1) = lumpedMTarget (2, 2) = lumpedMTarget (6, 6) = lumpedMTarget (7, 7) = lumpedMTarget (8, 8) = 2.7475E-07;

    Matrixd lumpedMTorsionTarget = lumpedMTest;
    lumpedMTorsionTarget (3, 3) = lumpedMTorsionTarget (9, 9) = 3.72057265E-07;

    const Matrixd coupledMTarget (12, 12, {
    2.28958333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.57916667E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,
    0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,
    0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,
    4.57916667E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.28958333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.87833333E-07,
    0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,
    0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07
    });

    const Matrixd consistentMTarget (12, 12, {
    1.83166667E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  9.15833333E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,
    0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.48038177E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.24019088E-07,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,
    0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,
    9.15833333E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.83166667E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.87833333E-07,
    0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.24019088E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.48038177E-07,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,
    0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07
    });

    // Check
    ASSERT_TRUE (lumpedMTarget * factor == lumpedMTest * factor);
    ASSERT_TRUE (lumpedMTorsionTarget * factor == lumpedMTorsionTest * factor);
    ASSERT_TRUE (coupledMTarget * factor == coupledMTest * factor);
    ASSERT_TRUE (consistentMTarget * factor == consistentMTest * factor);
}

// 单元坐标系下包含非结构质量的单元质量矩阵
TEST_F (BarCalculatorTest, ComputeLocalMassWithNSM)
{
    // Input data
    const Real L{10.0};

    IsotropicMaterial material;
    material.m_rho = 1.0;

    PBARData property;
    property.m_A   = 7.0;
    property.m_I1  = 7.145833;
    property.m_I2  = 2.333333;
    property.m_NSM = 1.0;
    bool hasTorsion{false};
    bool isBARMASS{false};

    // Test data
    auto lumpedMTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    hasTorsion = true;
    auto lumpedMTorsionTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    auto coupledMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    isBARMASS = true;
    auto consistentMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);

    // Target data
    Matrixd lumpedMTarget (12, 12);
    lumpedMTarget (0, 0) = lumpedMTarget (1, 1) = lumpedMTarget (2, 2) = lumpedMTarget (6, 6) = lumpedMTarget (7, 7) = lumpedMTarget (8, 8) = 40.0;

    Matrixd lumpedMTorsionTarget = lumpedMTest;
    lumpedMTorsionTarget (3, 3) = lumpedMTorsionTarget (9, 9) = 4.73958300E+01;

    const Matrixd coupledMTarget (12, 12, {
    3.33333333E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  6.66666667E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.47619048E+01,
    0.00000000E+00,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.47619048E+01,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00,
    0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -5.71428571E+01,
    6.66666667E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  3.33333333E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.19047619E+01,
    0.00000000E+00,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00, -2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,
    0.00000000E+00, -2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.61904762E+01
    });

    const Matrixd consistentMTarget (12, 12, {
    2.66666667E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.33333333E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.47619048E+01,
    0.00000000E+00,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  3.15972200E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.57986100E+01,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.47619048E+01,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00,
    0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -5.71428571E+01,
    1.33333333E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.66666667E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  1.02857143E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -4.19047619E+01,
    0.00000000E+00,  0.00000000E+00,  1.02857143E+01,  0.00000000E+00, -2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.97142857E+01,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.57986100E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  3.15972200E+01,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  2.47619048E+01,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  4.19047619E+01,  0.00000000E+00,  7.61904762E+01,  0.00000000E+00,
    0.00000000E+00, -2.47619048E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -5.71428571E+01,  0.00000000E+00, -4.19047619E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.61904762E+01
    });

    // Check
    ASSERT_TRUE (lumpedMTarget * factor == lumpedMTest * factor);
    ASSERT_TRUE (lumpedMTorsionTarget * factor == lumpedMTorsionTest * factor);
    ASSERT_TRUE (coupledMTarget * factor == coupledMTest * factor);
    ASSERT_TRUE (consistentMTarget * factor == consistentMTest * factor);
}

// 单元坐标系下单元耦合质量矩阵的自由度释放
TEST_F (BarCalculatorTest, ComputeLocalMassWithPin)
{
    // Input data
    const Real L{10.0};

    CBARData element;
    element.m_pb = 14;

    IsotropicMaterial material;
    material.m_E   = 210000.0;
    material.m_nu  = 0.3;
    material.m_G   = material.m_E / (2.0 * (1.0 + material.m_nu));
    material.m_rho = 7.85E-9;

    PBARData property;
    property.m_A   = 7.0;
    property.m_I1  = 7.145833;
    property.m_I2  = 2.333333;
    property.m_I12 = 0.0;
    property.m_J   = 6.003188;
    property.m_K1  = 0.833333;
    property.m_K2  = 0.833333;

    bool isBARMASS{false};

    // Test data
    auto coupledMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    isBARMASS = true;
    auto consistentMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    Matrixd kl (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kl);
    std::vector<std::size_t> freeDofSet;
    for (const auto& dof : element.GetDofSetOfPaRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof));
    }
    for (const auto& dof: element.GetDofSetOfPbRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof) + 6);
    }
    BarCalculator::ReleaseMassDof (freeDofSet, kl, _OUT coupledMTest);
    BarCalculator::ReleaseMassDof (freeDofSet, kl, _OUT consistentMTest);

    // Target data
    const Matrixd coupledMTarget (12, 12, {
    5.49500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,
    0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,
    0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.87833333E-07,
    0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,
    0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07
    });

    const Matrixd consistentMTarget (12, 12, {
    5.49500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,
    0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  7.44114531E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,
    0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  7.06500000E-08,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -2.87833333E-07,
    0.00000000E+00,  0.00000000E+00,  7.06500000E-08,  0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.04100000E-07,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    0.00000000E+00,  0.00000000E+00,  1.70083333E-07,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  2.87833333E-07,  0.00000000E+00,  5.23333333E-07,  0.00000000E+00,
    0.00000000E+00, -1.70083333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -3.92500000E-07,  0.00000000E+00, -2.87833333E-07,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  5.23333333E-07
    });

    // Check
    ASSERT_TRUE (coupledMTarget * factor == coupledMTest * factor);
    ASSERT_TRUE (consistentMTarget * factor == consistentMTest * factor);
}

TEST_F (BarCalculatorTest, ComputeGlobalMassWithPinOffsetNSM)
{
    // Input data
    const std::vector<NodeData> nodeSet{{6, 0, {0.0, 0.0, 0.0}, 0}, {8, 0, {3.0, 4.0, 12.0}, 0}};
    const NodeData orienNode;

    CBARData element;
    element.m_id   = 2;
    element.m_pId  = 2;
    element.m_g[0] = 6;
    element.m_g[1] = 8;
    element.m_flag = 65;  // BGG
    element.m_x1   = 3.0;
    element.m_x2   = 1.0;
    element.m_x3   = 4.0;
    element.m_pb   = 1;
    element.m_w1a  = 2.0;
    element.m_w2a  = 8.0;
    element.m_w3a  = 4.0;
    element.m_w1b  = 6.0;
    element.m_w2b  = 9.0;
    element.m_w3b  = 11.0;

    IsotropicMaterial material;
    material.m_E   = 210000.0;
    material.m_nu  = 0.3;
    material.m_G   = material.m_E / (2.0 * (1.0 + material.m_nu));
    material.m_rho = 1.0;

    PBARData property;
    property.m_A   = 7.0;
    property.m_I1  = 7.145833;
    property.m_I2  = 2.333333;
    property.m_I12 = 0.0;
    property.m_NSM = 1.0;
    property.m_J   = 6.003188;
    property.m_K1  = 0.833333;
    property.m_K2  = 0.833333;

    bool hasTorsion{false};
    bool isBARMASS{false};

    // Test data
    Real L;
    Matrixd T (3, 3);
    Matrixd TEle (12, 12);
    BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
    Matrixd kl (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kl);
    auto lumpedM414Test = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    auto lumpedMTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    hasTorsion = true;
    auto lumpedMTorsionTest = BarCalculator::ComputeLumpedMass (hasTorsion, L, material.m_rho, property);
    auto coupledMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    isBARMASS = true;
    auto consistentMTest = BarCalculator::ComputeCoupledMass (isBARMASS, L, material.m_rho, property);
    std::vector<std::size_t> freeDofSet;
    for (const auto& dof : element.GetDofSetOfPaRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof));
    }
    for (const auto& dof: element.GetDofSetOfPbRelease ())
    {
        freeDofSet.push_back (static_cast<std::size_t> (dof) + 6);
    }
    BarCalculator::ReleaseMassDof (freeDofSet, kl, _OUT coupledMTest);
    BarCalculator::ReleaseMassDof (freeDofSet, kl, _OUT consistentMTest);
    lumpedMTest  = TEle.Transpose () * lumpedMTest * TEle;
    lumpedMTorsionTest  = TEle.Transpose () * lumpedMTorsionTest * TEle;
    coupledMTest = TEle.Transpose () * coupledMTest * TEle;
    consistentMTest = TEle.Transpose () * consistentMTest * TEle;

    // Target data
    Matrixd lumpedM414Target (12, 12);
    lumpedM414Target (0, 0) = lumpedM414Target (1, 1) = lumpedM414Target (2, 2) = lumpedM414Target (6, 6) = lumpedM414Target (7, 7) = lumpedM414Target (8, 8) = 8.34266145E+01;

    const Matrixd lumpedMTarget (12, 12, {
     8.34266145E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  3.33706458E+02, -6.67412916E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  8.34266145E+01,  0.00000000E+00, -3.33706458E+02,  0.00000000E+00,  1.66853229E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  6.67412916E+02, -1.66853229E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00, -3.33706458E+02,  6.67412916E+02,  6.67412916E+03, -1.33482583E+03, -6.67412916E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     3.33706458E+02,  0.00000000E+00, -1.66853229E+02, -1.33482583E+03,  1.66853229E+03, -2.66965166E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    -6.67412916E+02,  1.66853229E+02,  0.00000000E+00, -6.67412916E+02, -2.66965166E+03,  5.67300978E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  9.17692759E+02, -7.50839530E+02,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  0.00000000E+00, -9.17692759E+02,  0.00000000E+00,  5.00559687E+02,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  7.50839530E+02, -5.00559687E+02,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -9.17692759E+02,  7.50839530E+02,  1.68521761E+04, -4.50503718E+03, -5.50615655E+03,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  9.17692759E+02,  0.00000000E+00, -5.00559687E+02, -4.50503718E+03,  1.30979785E+04, -8.25923483E+03,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -7.50839530E+02,  5.00559687E+02,  0.00000000E+00, -5.50615655E+03, -8.25923483E+03,  9.76091389E+03
    });

    const Matrixd lumpedMTorsionTarget (12, 12, {
     8.34266145E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  3.33706458E+02, -6.67412916E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  8.34266145E+01,  0.00000000E+00, -3.33706458E+02,  0.00000000E+00,  1.66853229E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  6.67412916E+02, -1.66853229E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00, -3.33706458E+02,  6.67412916E+02,  6.68526419E+03, -1.32687223E+03, -6.37189249E+02,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     3.33706458E+02,  0.00000000E+00, -1.66853229E+02, -1.32687223E+03,  1.67421343E+03, -2.64806333E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
    -6.67412916E+02,  1.66853229E+02,  0.00000000E+00, -6.37189249E+02, -2.64806333E+03,  5.75504545E+03,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  9.17692759E+02, -7.50839530E+02,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  0.00000000E+00, -9.17692759E+02,  0.00000000E+00,  5.00559687E+02,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  8.34266145E+01,  7.50839530E+02, -5.00559687E+02,  0.00000000E+00,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -9.17692759E+02,  7.50839530E+02,  1.68633112E+04, -4.49708358E+03, -5.47593289E+03,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  9.17692759E+02,  0.00000000E+00, -5.00559687E+02, -4.49708358E+03,  1.31036596E+04, -8.23764650E+03,
     0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00,  0.00000000E+00, -7.50839530E+02,  5.00559687E+02,  0.00000000E+00, -5.47593289E+03, -8.23764650E+03,  9.84294956E+03
    });

    const Matrixd coupledMTarget (12, 12, {
     7.37880322E+01,  8.43855411E+00,  3.20665056E+01,  2.22777828E+02,  3.97077807E+02, -6.17126805E+02,  1.90360630E+01, -1.72606789E+00, -6.55905796E+00, -4.00447749E+01,  1.50625451E+02, -1.55858450E+02,
     8.43855411E+00,  6.80015951E+01,  2.29046469E+01, -2.54827895E+02, -1.20550773E+01,  1.29674275E+02, -1.72606789E+00,  2.02196524E+01, -4.68504140E+00, -1.66455959E+02,  9.12350168E+00,  1.00700992E+02,
     3.20665056E+01,  2.29046469E+01,  1.49011715E+02,  1.14417478E+03, -2.30936924E+02, -2.10722751E+02, -6.55905796E+00, -4.68504140E+00,  3.64940067E+00,  5.85575380E+01, -5.78945087E+01,  3.09212733E+01,
     2.22777828E+02, -2.54827895E+02,  1.14417478E+03,  1.17999261E+04, -1.94229214E+03, -2.74794391E+03, -4.55681922E+01, -2.16484530E+02,  7.37578945E+01,  1.98603015E+03, -6.12871910E+02, -5.85676752E+02,
     3.97077807E+02, -1.20550773E+01, -2.30936924E+02, -1.94229214E+03,  3.48830100E+03, -3.52649368E+03,  1.87387957E+02,  2.46581126E+00, -6.96865661E+01, -5.60944665E+02,  1.52593990E+03, -1.45518482E+03,
    -6.17126805E+02,  1.29674275E+02, -2.10722751E+02, -2.74794391E+03, -3.52649368E+03,  5.78591106E+03, -1.81563163E+02,  9.03993808E+01,  4.31023809E+01, -2.51710014E+02, -1.35758253E+03,  1.80938765E+03,
     1.90360630E+01, -1.72606789E+00, -6.55905796E+00, -4.55681922E+01,  1.87387957E+02, -1.81563163E+02,  5.49930708E+01, -4.98641834E+00, -1.89483897E+01, -1.15684905E+02,  5.52555427E+02, -4.81156492E+02,
    -1.72606789E+00,  2.02196524E+01, -4.68504140E+00, -2.16484530E+02,  2.46581126E+00,  9.03993808E+01, -4.98641834E+00,  5.84123291E+01, -1.35345641E+01, -5.98288007E+02,  2.63567826E+01,  3.34172222E+02,
    -6.55905796E+00, -4.68504140E+00,  3.64940067E+00,  7.37578945E+01, -6.96865661E+01,  4.31023809E+01, -1.89483897E+01, -1.35345641E+01,  1.05427131E+01,  2.00064967E+02, -2.10509047E+02,  8.93281228E+01,
    -4.00447749E+01, -1.66455959E+02,  5.85575380E+01,  1.98603015E+03, -5.60944665E+02, -2.51710014E+02, -1.15684905E+02, -5.98288007E+02,  2.00064967E+02,  6.77519454E+03, -1.97792585E+03, -2.08693662E+03,
     1.50625451E+02,  9.12350168E+00, -5.78945087E+01, -6.12871910E+02,  1.52593990E+03, -1.35758253E+03,  5.52555427E+02,  2.63567826E+01, -2.10509047E+02, -1.97792585E+03,  5.79896342E+03, -4.48512438E+03,
    -1.55858450E+02,  1.00700992E+02,  3.09212733E+01, -5.85676752E+02, -1.45518482E+03,  1.80938765E+03, -4.81156492E+02,  3.34172222E+02,  8.93281228E+01, -2.08693662E+03, -4.48512438E+03,  5.69265956E+03
    });

    const Matrixd consistentMTarget (12, 12, {
     7.37880322E+01,  8.43855411E+00,  3.20665056E+01,  2.22777828E+02,  3.97077807E+02, -6.17126805E+02,  1.90360630E+01, -1.72606789E+00, -6.55905796E+00, -4.00447749E+01,  1.50625451E+02, -1.55858450E+02,
     8.43855411E+00,  6.80015951E+01,  2.29046469E+01, -2.54827895E+02, -1.20550773E+01,  1.29674275E+02, -1.72606789E+00,  2.02196524E+01, -4.68504140E+00, -1.66455959E+02,  9.12350168E+00,  1.00700992E+02,
     3.20665056E+01,  2.29046469E+01,  1.49011715E+02,  1.14417478E+03, -2.30936924E+02, -2.10722751E+02, -6.55905796E+00, -4.68504140E+00,  3.64940067E+00,  5.85575380E+01, -5.78945087E+01,  3.09212733E+01,
     2.22777828E+02, -2.54827895E+02,  1.14417478E+03,  1.18073495E+04, -1.93698974E+03, -2.72779480E+03, -4.55681922E+01, -2.16484530E+02,  7.37578945E+01,  1.98974183E+03, -6.10220711E+02, -5.75602196E+02,
     3.97077807E+02, -1.20550773E+01, -2.30936924E+02, -1.93698974E+03,  3.49208843E+03, -3.51210146E+03,  1.87387957E+02,  2.46581126E+00, -6.96865661E+01, -5.58293466E+02,  1.52783362E+03, -1.44798871E+03,
    -6.17126805E+02,  1.29674275E+02, -2.10722751E+02, -2.72779480E+03, -3.51210146E+03,  5.84060151E+03, -1.81563163E+02,  9.03993808E+01,  4.31023809E+01, -2.41635458E+02, -1.35038642E+03,  1.83673287E+03,
     1.90360630E+01, -1.72606789E+00, -6.55905796E+00, -4.55681922E+01,  1.87387957E+02, -1.81563163E+02,  5.49930708E+01, -4.98641834E+00, -1.89483897E+01, -1.15684905E+02,  5.52555427E+02, -4.81156492E+02,
    -1.72606789E+00,  2.02196524E+01, -4.68504140E+00, -2.16484530E+02,  2.46581126E+00,  9.03993808E+01, -4.98641834E+00,  5.84123291E+01, -1.35345641E+01, -5.98288007E+02,  2.63567826E+01,  3.34172222E+02,
    -6.55905796E+00, -4.68504140E+00,  3.64940067E+00,  7.37578945E+01, -6.96865661E+01,  4.31023809E+01, -1.89483897E+01, -1.35345641E+01,  1.05427131E+01,  2.00064967E+02, -2.10509047E+02,  8.93281228E+01,
    -4.00447749E+01, -1.66455959E+02,  5.85575380E+01,  1.98974183E+03, -5.58293466E+02, -2.41635458E+02, -1.15684905E+02, -5.98288007E+02,  2.00064967E+02,  6.78261790E+03, -1.97262345E+03, -2.06678751E+03,
     1.50625451E+02,  9.12350168E+00, -5.78945087E+01, -6.10220711E+02,  1.52783362E+03, -1.35038642E+03,  5.52555427E+02,  2.63567826E+01, -2.10509047E+02, -1.97262345E+03,  5.80275085E+03, -4.47073216E+03,
    -1.55858450E+02,  1.00700992E+02,  3.09212733E+01, -5.75602196E+02, -1.44798871E+03,  1.83673287E+03, -4.81156492E+02,  3.34172222E+02,  8.93281228E+01, -2.06678751E+03, -4.47073216E+03,  5.74735000E+03
    });

    // Check
    ASSERT_TRUE (lumpedM414Target == lumpedM414Test);
    ASSERT_TRUE (lumpedMTarget == lumpedMTest);
    ASSERT_TRUE (lumpedMTorsionTarget == lumpedMTorsionTest);
    ASSERT_TRUE (coupledMTarget == coupledMTest);
    ASSERT_TRUE (consistentMTarget == consistentMTest);
}

TEST_F (BarCalculatorTest, ComputeGlobalElementInternalForceStrainStress)
{
    // Input data
    const std::vector<NodeData> nodeSet{{7, 0, {0.0, 0.0, 0.0}, 0}, {8, 0, {100.0, 0.0, 0.0}, 0}};
    const NodeData orienNode;

    CBARData element;
    element.m_g[0] = 7;
    element.m_g[1] = 8;
    element.m_flag = 1;  // GGG
    element.m_x1   = 3.0;
    element.m_x2   = 2.0;
    element.m_x3   = 1.0;

    IsotropicMaterial material;
    material.m_id = 2;
    material.m_E  = 210000.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2.0 * (1.0 + material.m_nu));

    const PBARData property{
        2,        2,        11.0,     28.95644, 8.303030, 14.60417, 0.0,  0.0,
        1.75,     1.0,      -1.75,    1.0,      -1.75,    -1.0,     1.75, -1.0,
        0.318182, 0.454545, -7.15909, 0};

    const Matrixd eleDisplacement (12, 1, {
    0.0,                  0.0,                  0.0,               0.0,                   0.0,                  0.0,
    4.329004329004329E-5, 0.015898800082727524, 0.689996925379203, 8.477683004889961E-5, -0.010287280973703707, 2.544100563841669E-4
    });

    // Test data
    Real L;
    Matrixd T (3, 3);
    Matrixd TEle (12, 12);
    BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
    Matrixd kl (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kl);
    Matrixd eleInternalForceTest (12, 1);
    BarCalculator::ComputeElementInternalForce (kl, TEle, eleDisplacement, _OUT eleInternalForceTest);
    Matrixd eleStressTest (9, 1);
    BarCalculator::ComputeElementStress (eleInternalForceTest, property, _OUT eleStressTest);
    Matrixd eleStrainTest (9, 1);
    BarCalculator::ComputeElementStrain (eleStressTest, material.m_E, _OUT eleStrainTest);

    // Target data
    const Matrixd eleInternalForceTarget (12, 1, {
    -1.000000E+00, -3.130495E+00, -1.788854E+00, -1.000000E+00,  1.753077E+02, -3.134967E+02,
     1.000000E+00,  3.130495E+00,  1.788854E+00,  1.000000E+00,  3.577709E+00,  4.472136E-01
    });
    const Matrixd eleStressTarget (9, 1, {
    9.090909E-02, -7.438758E+01, -3.008401E+00,  7.438758E+01,  3.008401E+00,
                   7.332994E-01,  3.281177E-01, -7.332994E-01, -3.281177E-01
    });
    const Matrixd eleStrainTarget (9, 1, {
    4.329004E-07, -3.542266E-04, -1.432572E-05,  3.542266E-04,  1.432572E-05,
                   3.491902E-06,  1.562465E-06, -3.491902E-06, -1.562465E-06
    });

    // Check
    ASSERT_TRUE (eleInternalForceTarget == eleInternalForceTest);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
    ASSERT_TRUE (eleStrainTarget == eleStrainTest);
}

TEST_F (BarCalculatorTest, ComputeGlobalElementInternalForceStrainStressWithPin)
{
    // Input data
    const std::vector<NodeData> nodeSet{{7, 0, {0.0, 0.0, 0.0}, 0}, {8, 0, {100.0, 0.0, 0.0}, 0}};
    const NodeData orienNode;

    CBARData element;
    element.m_g[0] = 7;
    element.m_g[1] = 8;
    element.m_flag = 1;  // GGG
    element.m_pb   = 1;
    element.m_x1   = 3.0;
    element.m_x2   = 2.0;
    element.m_x3   = 1.0;

    IsotropicMaterial material;
    material.m_id = 2;
    material.m_E  = 210000.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2.0 * (1.0 + material.m_nu));

    const PBARData property{
        2,        2,        11.0,     28.95644, 8.303030, 14.60417, 0.0,  0.0,
        1.75,     1.0,      -1.75,    1.0,      -1.75,    -1.0,     1.75, -1.0,
        0.318182, 0.454545, -7.15909, 0};

    const Matrixd eleDisplacement (12, 1, {
    0.0,   0.0,                0.0,               0.0,                   0.0,                  0.0,
    0.0, 0.015898800082727524, 0.689996925379203, 8.477683004889961E-5, -0.010287280973703707, 2.544100563841669E-4
    });

    // Test data
    Real L;
    Matrixd T (3, 3);
    Matrixd TEle (12, 12);
    BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
    Matrixd kl (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kl);
    Matrixd eleInternalForceTest (12, 1);
    BarCalculator::ComputeElementInternalForce (kl, TEle, eleDisplacement, _OUT eleInternalForceTest);
    Matrixd eleStressTest (9, 1);
    BarCalculator::ComputeElementStress (eleInternalForceTest, property, _OUT eleStressTest);
    Matrixd eleStrainTest (9, 1);
    BarCalculator::ComputeElementStrain (eleStressTest, material.m_E, _OUT eleStrainTest);

    // Target data
    const Matrixd eleInternalForceTarget (12, 1, {
    0.0, -3.130495E+00, -1.788854E+00, -1.000000E+00,  1.753077E+02, -3.134967E+02,
    0.0,  3.130495E+00,  1.788854E+00,  1.000000E+00,  3.577709E+00,  4.472136E-01
    });
    const Matrixd eleStressTarget (9, 1, {
    0.0, -7.438758E+01, -3.008401E+00,  7.438758E+01,  3.008401E+00,
          7.332994E-01,  3.281177E-01, -7.332994E-01, -3.281177E-01
    });
    const Matrixd eleStrainTarget (9, 1, {
    0.0, -3.542266E-04, -1.432572E-05,  3.542266E-04,  1.432572E-05,
          3.491902E-06,  1.562465E-06, -3.491902E-06, -1.562465E-06
    });

    // Check
    ASSERT_TRUE (eleInternalForceTarget == eleInternalForceTest);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
    ASSERT_TRUE (eleStrainTarget == eleStrainTest);
}

TEST_F (BarCalculatorTest, ComputeGlobalElementInternalForceStrainStressWithOffset)
{
    // Input data
    const std::vector<NodeData> nodeSet{{7, 0, {0.0, 0.0, 0.0}, 0}, {8, 0, {100.0, 0.0, 0.0}, 0}};
    const NodeData orienNode;

    CBARData element;
    element.m_g[0] = 7;
    element.m_g[1] = 8;
    element.m_flag = 1;  // GGG
    element.m_x1   = 3.0;
    element.m_x2   = 2.0;
    element.m_x3   = 1.0;
    element.m_w1a  = 3.0;
    element.m_w2a  = 4.0;
    element.m_w3a  = 5.0;
    element.m_w1b  = 4.0;
    element.m_w2b  = 5.0;
    element.m_w3b  = 3.0;

    IsotropicMaterial material;
    material.m_id = 2;
    material.m_E  = 210000.0;
    material.m_nu = 0.3;
    material.m_G  = material.m_E / (2.0 * (1.0 + material.m_nu));

    const PBARData property{
        2,        2,        11.0,     28.95644, 8.303030, 14.60417, 0.0,  0.0,
        1.75,     1.0,      -1.75,    1.0,      -1.75,    -1.0,     1.75, -1.0,
        0.318182, 0.454545, -7.15909, 0};

    const Matrixd eleDisplacement (12, 1, {
    0.0,                 0.0,                   0.0,                 0.0,                  0.0,                  0.0,
    0.04262842065574956, 0.0021159649366413492, 0.6289023125003526, -5.75690854113826E-4, -0.009646961047349712, 1.0542666016150824E-4
    });

    // Test data
    Real L;
    Matrixd T (3, 3);
    Matrixd TEle (12, 12);
    BarCalculator::ComputeElementTrans (element, nodeSet, orienNode, _OUT L, _OUT T, _OUT TEle);
    Matrixd kl (12, 12);
    BarCalculator::ComputeStiffness (L, material, property, _OUT kl);
    Matrixd eleInternalForceTest (12, 1);
    BarCalculator::ComputeElementInternalForce (kl, TEle, eleDisplacement, _OUT eleInternalForceTest);
    Matrixd eleStressTest (9, 1);
    BarCalculator::ComputeElementStress (eleInternalForceTest, property, _OUT eleStressTest);
    Matrixd eleStrainTest (9, 1);
    BarCalculator::ComputeElementStrain (eleStressTest, material.m_E, _OUT eleStrainTest);

    // Target data
    const Matrixd eleInternalForceTarget (12, 1, {
    -9.601608E-01, -3.182858E+00, -1.716831E+00,  7.859460E+00,  1.633521E+02, -3.148084E+02,
     9.601608E-01,  3.182858E+00,  1.716831E+00, -7.859460E+00,  1.009033E+01, -6.739004E+00
    });
    const Matrixd eleStressTarget (9, 1, {
    8.728734E-02, -7.191615E+01, -1.919051E+00,  7.191615E+01,  1.919051E+00,
                   2.985409E+00,  6.136669E-01, -2.985409E+00, -6.136669E-01
    });
    const Matrixd eleStrainTarget (9, 1, {
    4.156540E-07, -3.424579E-04, -9.138340E-06,  3.424579E-04,  9.138340E-06,
                   1.421623E-05,  2.922224E-06, -1.421623E-05, -2.922224E-06
    });

    // Check
    ASSERT_TRUE (eleInternalForceTarget == eleInternalForceTest);
    ASSERT_TRUE (eleStressTarget == eleStressTest);
    ASSERT_TRUE (eleStrainTarget == eleStrainTest);
}